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Abstract 

The relatively large measured value of ^13 has opened up the possibility of determining the neutrino 
mass hierarchy through earth matter effects. Amongst the current accelerator-based experiments 
only NOvA has a long enough baseline to observe earth matter effects. However, NOvA is plagued 
with uncertainty on the knowledge of the true value of 6c p, and this could drastically reduce its 
sensitivity to the neutrino mass hierarchy. The earth matter effect on atmospheric neutrinos on 
the other hand is almost independent of Sep- The 50 kton magnetized Iron CALorimeter at the 
India-based Neutrino Observatory (ICALQINO) will be observing atmospheric neutrinos. The charge 
identification capability of this detector gives it an edge over others for mass hierarchy determination 
through observation of earth matter effects. We study in detail the neutrino mass hierarchy sensitivity 
of the data from this experiment simulated using the Nuance based generator developed for ICAL@INO 
and folded with the detector resolutions and efficiencies obtained by the INO collaboration from a full 
Geant4-based detector simulation. The data from ICAL@INO is then combined with simulated data 
from T2K, NOvA, Double Chooz, RENO and Daya Bay experiments and a combined sensitivity study 
to the mass hierarchy is performed. With 10 years of ICAL@INO data combined with T2K, NOvA 
and reactor data, one could get about 2.2a — 5.5(7 discovery of the neutrino mass hierarchy, depending 
on the true value of sin^ 6*23 [0.4 - 0.6], sin^ 26*13 [0.08 - 0.12] and 8cp [0 - 2-k\. 



*email:aiiushree@hri . res . in 
^email:tarak(§tif r . res . in 
■''email: sandhyaOhri . res . in 

1 



1 Introduction 



The past year has been very exciting in the field of neutrino physics. Results from as many as five 
different experiments finally confirmed that the mixing angle ^13 is non-zero and indeed Zarg^. The first 
hints came in 2011 from the T2K experiment [1], which with just 6 events collected before its forced 
temporary halt due to the great Japan earthquake, inferred a value of sin^ 2^13 ~ 0.1, with a zero value 
excluded at the 2.5a C.L.. This was corroborated with similar hints from the MINOS [2J and the Double 
Chooz [3] experiments. Finally, the Daya Bay experiment put all doubts to rest in March 2012, when 
with just 55 live days of their data they established a non-zero ^13 at the 5.2a C.L. with the best-fit 
sin^ 26113 = 0.092 ± 0.016(stat) ± 0.005(syst) g]. These results were further strengthened by the RENO 
experiment which provided an independent confirmation of ^13 7^ at the 4.9cr C.L. in April 2012, with a 
best-fit sin^ 2^13 = 0.113ib0.013(stat) ib0.019(syst) [5j. All these experiments confirmed their first results 
with more statistics at the Neutrino Conference in Kyoto, Japan [6J. A global analysis of the world 
neutrino data now provides a more than lOcr signal for a non-zero ^13, with best-fit sin^ 2^13 ~ 0.1 it 0.01 

mum. 

The relatively large value of ^13 has opened up the possibility of answering the two other major elusive 
issues in neutrino oscillation physics - CP violation in the lepton sector, and the sign of Am'^i, aka, the 
neutrino mass hierarch}d. In order to address these two issues a huge effort worldwide is underway and 
next generation experiments are being proposed. In the light of the large measured value of sin^ 2^13, one 
needs to take another look at the physics reach of these prospective experimental endeavors to allow one 
to make prudent choices for the design of the next generation long baseline experiments. While the large 
value of sin^ 2^13 would necessarily improve the prospects of determining the neutrino mass hierarchy 
through observations of earth matter effects in neutrino beam experiments, the same is not necessarily 
true for the measurement of CP violation [10]. The latter effect appears through the sub-dominant term 
in the i/^ — )• oscillation probability, and hence is much more subtle. Therefore, it is pertinent to ask 
at this point if one could measure the neutrino mass hierarchy elsewhere [11] . One could then optimize 
the long baseline experiment mainly to measure CP violation. In particular, with sin^ 2^13 ~ 0.1 one 
could use the earth matter effects in atmospheric neutrinos in order to look for, and possibly identify the 
neutrino mass hierarchy. 

Prospects of detecting earth matter effects in atmospheric neutrinos have been discussed in detail in 
the past by many authors [12|-[31|. In particular, the possibility of measuring the neutrino mass hierarchy 
in magnetized iron calorimeters has been considered before. In this paper, we expound the reach of 
atmospheric neutrino measurements at the magnetized Iron CALorimeter (ICAL) detector to be built at 
the India-based Neutrino Observatory (INO) for the determination of the mass hierarchy. Henceforth we 
will designate this experiment as ICAL@INO. The main new elements in this work are the simulation 
of atmospheric neutrino events in ICAL@INO and its detailed physics analysis, for which we use the 
tools developed by the INO collaboration for this experiment. Most past works on the mass hierarchy 
determination with magnetized iron calorimeters performed the analysis in terms of the neutrino energy 
and zenith angle bins, with some assumed fixed values for the resolutions and efficiencies. In this work, 
we perform our analysis in terms of the muons, which are binned in reconstructed energy and zenith angle 
bins. For obtaining the reconstructed energy and zenith angle of the muons we use, for the first time, 
the detailed efficiencies and resolution functions obtained by the INO collaboration from the simulation 

^ "Large" here implies that the measured value of ^13 turned out to be larger than expected. In fact, it was found to lie 
just below the earlier upper bound on this parameter. However, it is still smaller compared to the other two mixing angles 
023 and 9i2. 

^We define mass squared differences as Am^j — m| — . 
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codes developed for ICAL. In particular, we use the Nuance V3.000 [32] based atmospheric neutrino 
event generator used by the INO collaboration to simulate the events in the detector in the absence of 
oscillations. This raw Nuance data is then folded with the oscillation code via a reweighting algorithm 
and the oscillated event spectrum generated. These events are then folded with the resolutions and 
efficiencies obtained by the INO collaboration from a full Geant4-based detector Monte Carlo developed 
to simulate the IC AL@INO detector [33] . Separate papers will appear describing in detail the simulation 
codes and the results from these studies performed by the collaboration [33\ I34j . Here we use these 
results to do the physics analysis of simulated atmospheric neutrino events in ICAL@INO. In particular, 
we study the reach of the experiment in pinning down the true neutrino mass hierarchy as a function of 
the number of years of running of the experiment. We quantitatively show how this sensitivity depends 
on the uncertainties in the measurement of the other oscillation parameters, especially lAm^^l, 623 and 
^13. Since all three of these are expected to be measured to a remarkable precision by the current reactor 
and accelerator-based long baseline experiments, we include in our analysis the simulated data from the 
full run of all of these experiments and show the joint sensitivity reach to the neutrino mass hierarchy. 
We find that while these reactor and accelerator-based neutrino oscillation experiments themselves have 
very limited sensitivity to the neutrino mass hierarchy, they still have a crucial role to play in this effort, 
since they constrain lAmgj^l, 623 and ^13, which in turn improves the statistical significance with which 
ICAL@INO can determine the neutrino mass hierarchy. We also show that the reach of ICAL@INO 
for the neutrino mass hierarchy is nearly independent of Sep, which at the moment is totally unknown 
and which is also the hardest amongst all oscillation parameters to be measured. The long baseline 
experiments, on the other hand are bogged down by the uncertainty in the true value of 5c p, making it 
difficult to measure the neutrino mass hierarchy from these experiments alone [35]. We will show that 
ICAL@INO will provide a remarkable complementarity in this direction. 

The paper is organized as follows. In section 2 we briefiy review the current status of the already 
measured neutrino oscillation parameters and give their benchmark values at which we simulate the 
projected data used in our analysis. We give the ranges of these parameters over which they are allowed 
to vary in our statistical fits. In section 3 we describe the experiments ICAL@INO, Double Chooz, Daya 
Bay, RENO, T2K, and NOi^A and spell out the experimental specifications used in our analysis for each 
one of them. In section 4 we present the simulated events at ICAL@INO and show the impact of the 
efficiencies and resolutions for the muons in ICAL obtained by the INO collaboration. In section 5 we 
give the details of our statistical analysis tool. In section 6 we present our main results. The impact of 
systematic uncertainties on the mass hierarchy sensitivity of ICAL@INO is discussed in section 7 and 
that of the true value of 623 is studied in section 8. In section 9 we explore the effect of 6c p value on the 
mass hierarchy sensitivity in ICAL@INO and NOz/A. We end with our conclusions in section 10. 

2 Neutrino Oscillation Paramaters 

We start with a brief overview of the current status of the neutrino oscillation parameters [3 [H [9j and 
their best-fit values which we use for simulating events in the various experiments. The solar neutrino 
parameters Am^i and sin^ 612 are now determined to extremely good precision by the joint analysis of 
the solar |36j and KamLAND [37] data. Observation of matter effects in solar neutrinos has also nailed 
down the sign of Am^i to be positive. The current best-fit from a global analysis [8j is 

Am^i = (7.54 ± 0.26) x W'^eV'^, sin^ 6*12 = 0.31 ± 0.02 . (1) 
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Parameter 


True value used in data 


3a range used in fit 




7.5 X 10"^ eV^ 


[7.0 - 8.0] X 10-5 eV^ 


sin^ 6*12 


0.3 


[0.265 - 0.33] 




2.4 X 10^3 eV^ 


[2.1 - 2.6] X lO^^ eV^ 


Sep 





[0 - 2tt] 


sin^ 6*23 
sin^ 26^13 


0.4, 0.5, 0.6 
0.08, 0.1, 0.12 


sin^ 6^23 (true) ±0.1 
sin^ 26li3(true)±0.03 



Table 1: Benchmark true values of oscillation parameters used in the simulations, unless otherwise stated. 
The range over which they are allowed to vary freely in the fit is also shown in the last column. For 
sin^ ^23 (true) and sin^ 2^13 (true) we use three benchmark values for simulating the data. 



Among the atmospheric neutrino parameters, sin^ 2^23 is mainly constrained by the Super-Kamiokande 
atmospheric neutrino data [38], while [Am^^j is predominantly determined by the MINOS long baseline 
disappearance data [39]. For [Aml^l the current best-fit is close to [Till] 

[Am^il = 2.47 X IQ-^eV^ . (2) 

With neutrino physics entering the precision era, it has become very important to define what is meant by 
the atmospheric neutrino mass splitting when one is doing a three- generation fit. The subtlety involved is 
the following. The value of the best-fit for the atmospheric mass squared difference depends on the mass 
hierarchy and definition used. In particular, it could be rather misleading to use an inconsistent definition 
for this parameter when doing mass hierarchy studies. For instance, if Am^^ is called the atmospheric 
neutrino mass squared difference and Am^^ > defined as normal hierarchy, then the absolute value of 
Am32 changes when one changes the hierarchy from normal {Am^i > 0) to inverted [Am'^i < 0). Since 
the three generation oscillation probability is sensitive to all oscillation frequencies, in this case one gets 
a rather large difference in the survival probability Pu^u^ between normal and inverted hierarchies even 
when ^13 is taken as zero and there are no ^13 driven earth matter effects. One needs to perform a careful 
marginalization over Am^i in this case to get rid of the spurious difference in Py^^u^ coming from this 
effect [ID]. Therefore, it is important to use a consistent definition for the mass squared difference in the 
analysis, especially in studies pertaining to observations of earth matter effects. In our study we use as 
the atmospheric mass squared difference, the quantity defined as [41j 

Am^fj = Am|]^ — (cos^ 612 — cos 5cp sin ^13 sin 2^12 tan ^23) Am^i , (3) 

where the other parameters are defined according to the convention used by the PDG. The normal 
hierarchy is then defined as Am^g > and the inverted hierarchy as Am^g < 0. Defining the mass 
squared difference by Eq. ([3|) is particularly convenient for mass hierarchy studies involving the probability 
Pu^Uf^ since it is almost same for Amg^ > (normal hierarchy) and Am^g < (inverted hierarchy) for 
^13 = 0. This ensures that there is no spurious contribution to the mass hierarchy sensitivity coming from 
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the difference between the oscihation frequencies for the normal and inverted hierarchies in experiments 
predominantly sensitive to the oscillation channel Pj/^jy^- However for this definition of the mass hierarchy, 
there is a difference between the frequencies involved in normal and inverted ordering for the other 
oscillation channels and marginalizing over Am^g then becomes very important for them. In our analysis 
we have paid special attention to the marginalizing procedure and have checked our mass hierarchy 
sensitivity to the definition used for the atmospheric mass squared difference and the mass hierarchy. We 
will quantify this issue later. 

The issue regarding the value of sin^ 623 and its correct octant is not yet settled. The Super- 
Kamiokande collaboration get the best-fit to their atmospheric zenith angle data at sin^ 2^23 = 0.98 
|38j . The MINOS collaboration best-fit sin^ 2^23 = 0.96 [39], where they included in their analysis the 
full MINOS data with 10.71 x lO^o POT for u^-heam, 3.36 x 10^° POT for P^-beams, as well as their 37.9 
kton-years data from atmospheric neutrinos. It is worth pointing out here that the Super-Kamiokande 
best-fit is close to maximal and even the MINOS collaboration results allow maximal mixing at la C.L. 
|39j . However, results from global analyses performed by groups outside the experimental collaborations 
have now started to show deviation from maximal mixing and prefer sin^ ^23 in the first octant [H [7] . 
These results though should be taken only as a possible "hint" as they require further investigation by 
the analysis using the full detector Monte Carlo of the experimental collaborations. 

Our analysis in this paper uses the full three generation oscillation probabilities without any ap- 
proximations. The data are simulated at a particular set of benchmark values chosen for the oscillation 
parameters, which we call "true value". These are summarized in Table [TJ The true values of Am^]^, 
sin^ 612 and lAm^gl are kept fixed throughout the paper. These quantities are now fairly well determined 
by the current global neutrino data and we choose our benchmark true values for these parameters to be 
close to their current best-fits, as discussed above. We will show results for a range of plausible values 
of sin^ 023(true) and sin^ 2^13 (true) since the earth matter effects are fairly sensitive to these parameters. 
While we have absolutely no knowledge on the value of 5(7p(true), the sensitivity of ICAL@INO does not 
depend much on the true value of this parameter. Therefore, (5cp(true) is also kept fixed at zero, unless 
otherwise stated. Only in section 8, where we study the impact of 6c p (true) on the mass hierarchy reach 
of the NOz/A experiment, will we show results as a function of the 6cp(tTue). 

In our fit, we allow the oscillation parameters to vary freely within their current 3a limits and the 
is minimized (marginalized) over them. The range over which these parameters are varied in the fit is 
shown in Table [TJ Since the ICAL@INO sensitivity does not depend much on Am^]^, sin^ ^12 and 6cp, 
we keep these fixed at their true values in the fit for the analysis of the ICAL@INO data. However, the 
ICAL@INO sensitivity to the mass hierarchy does depend on |Amgg|, sin^ ^23 and sin^ 2^13 and hence 
the xfno is marginalized over the 3a ranges of these parameters. The combined for the accelerator 
and reactor experiments depends on all the oscillation parameters, and so for them we marginalize the 

over the current 3a range of all the oscillation parameters. Since the range of |Amg£j|, sin^ ^23 and 
sin^ 2^13 will be severely constrained by the future accelerator and reactor data themselves, the best-fit 
for these in our global fits are mostly close to the true value taken for the data. However, none of the 
data sets included in our analysis has the potential to constrain Aml]^ and sin^0i2. Therefore, in order 
to take into account the fact that not all values of Am^x and sin^ 9i2 within their current 3a range are 
allowed with equal probability by the solar and KamLAND data, we impose a "prior" according to the 
following definition: 

^prior ^ y ( ' ) ' (^) 



5 



where the parameter pi is Am^i and sin^ 6*12 for i = 1 and 2 respectively, with o"Am|i ~ "^"^ ''^sin^ 6»i2 ~ 
4%. This Xprior added to the sum of the obtained from the analysis of the ICAL@INO simulated 
data and the prospective accelerator and reactor data, and the combined is then marginalized over all 
oscillation parameters. 

3 Neutrino Oscillation Experiments 

We give below a very brief description of the experiments whose simulated data wc use in this analysis 
to pin down the neutrino mass hierarchy. The ICAL@INO experiment is of course the focus of this work. 
We begin with a short overview of this experimental endeavor. We then move on to mention just the 
key features of the accelerator-based experiments T2K and NOi/A and the reactor experiments Double 
Chooz, Daya Bay and RENO. 

3.1 ICAL@ India- based Neutrino Observatory 

ICAL will be a 50 kton magnetized iron calorimeter at the INO laboratory in India and will soon go 
into construction in the Theni district in Southern India. It will be solid rectangular in shape with 
dimensions 48.4 m in length, 16 m in width and 14.4 m in height and will consist of three identical 
modules. The detector will have a layered structure with 150 layers of 5.6 cm iron slabs interleaved 
with glass Resistive Plate Chambers (RPC) acting as the active detector element. Each glass RFC to be 
deployed in ICAL@INO will be about 2 m x 2 m in size made up of two parallel glass electrodes separated 
by spacers to create a gap which is filled with a gas mixture of tetrafiouroethane and isobutane. This 
particular combination of gases chosen by the INO collaboration enables the RPC to be used in the 
avalanche mode wherein the arrival of a charged particle results in a Townsend avalanche through the 
gas volume. A total of about 28,000 units of such RPCs will be required for the complete detector. Eight 
such units will be arranged next to each other to form 16 m x 2 m road and each module will have eight 
such roads per layer. This signal in the RPC will be read by 3 cm wide pickup strips that will be laid 
orthogonal to each other (X and Y strips). The signal will then go to front end ASICs located at the end 
of the strips. There will be 64 strips long x direction and 64 strips along the y direction per RPC. Thus 
one needs about 3.7 million electronic readout channels for the full detector. The orthogonal X-Y readout 
strips form a grid such that muons traveling in the detector will trigger the RPCs in a particular 3 cm 
X 3 cm block which serves as a hit point. Since the muon will cross a number of RPCs, the hit points 
can be joined to reconstruct the long well defined muon track. In addition, being made of iron, the ICAL 
detector will be magnetized with a magnetic field strength of about 1.3 Tesla. Therefore as it travels, the 
IJ,~ bends in a direction opposite to that of the This gives ICAL an edge over other detectors since 
the magnetic field allows charge discrimination allowing it to distinguish between muon neutrinos and 
antineutrinos. On the other hand, ICAL is expected to not have very good sensitivity to electrons since 
the electron showers are mostly absorbed in the dense iron material. However, hadron showers can be seen 
in the detector and their energy measured. This makes this detector a calorimeter wherein the energy and 
momentum of the incoming (anti)neutrino can be reconstructed by adding the energy and momentum of 
the resultant muon and hadron(s). Therefore the energy and momentum resolution of ICAL is expected 
to be better than that of detectors which are insensitive to the energy and momentum of the final state 
hadrons. At least four types of large detectors for atmospheric neutrinos have been envisaged. While 
magnetized iron calorimeters such as ICAL@INO have excellent charge identification capabilities and 
good energy and angle resolution, they suflFer from diflaculty in observing electrons and have a higher 
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energy threshold for muons. The water Cherenkov detectors hke Hyper-Kamiokande fl2| observe both 
muons and electrons with very low energy threshold, but cannot be magnetized. Liquid argon detectors 
|301l43j have extremely good detector response but magnetization could still be a challenge for them |44j . 
The multi-megaton ice detector PINGU (IceCube extension for low-energy neutrinos) [35] could use its 
huge statistics to overcome its other drawbacks to return good sensitivity to the mass hierarchy |46j- The 
simulation of events in ICAL will be discussed in detail in the next section. 

3.2 Current Reactor and Accelerator Experiments 

For simulation of the current reactor and accelerator-based experiments we use the GLoBES software |47j . 
We have closely followed [35] for the analysis of the future accelerator and reactor data. The experiments 
that we include in our study are the following: 

• Double Chooz: The Double Chooz reactor experiment |48j has a liquid scintillator detector with 
fiducial mass of 8.3 tons placed at a distance of 1 km and 1.1 km from the two reactor cores of the 
Chooz reactor power plant, each with 4.27 GWth thermal power. Double Chooz has been taking 
data with just this far detector and have observed a positive signal for ^13 at 3.1cr C.L. O [6]. In 
addition to the far detector, this experiment will also have a near detector which will be identical 
to the far detector and will be placed at a distance of 470 m and 350 m respectively from the two 
reactor cores. Following [35j, we perform the analysis for an exposure of 3 years with both near and 
far detectors fully operational and with detector efficiency of 80% and reactor load factor of 78%. 
An uncorrelated systematic uncertainty of 0.6% is assumed. 

• RENO: The RENO antineutrino experiment [l9] is powered by the Yonggwang reactor plant in 
South Korea, with a total reactor power of 16.4 GWth, making it currently the most powerful reactor 
in the world behind the Kashiwazaki-Kariwa power plant in Japan (which is currently shutdown). 
This reactor complex has six reactor cores with first two having power 2.6 GWth while the last four 
with power 2.8 GWth, respectively. These reactor cores are arranged along a 1.5 km straight line 
separated from each other by equal distances. The near detector of this experiment has a fiducial 
mass of 15 tons and is situated at a distance of 669 m, 453 m, 307 m, 338 m, 515 m and 74 m 
respectively from the reactor cores. The far detector also has a fiducial mass of 15 tons and is placed 
in the opposite direction at a distance of 1.557 km, 1.457 km, 1.396 km, 1.382 km, 1.414 km and 
1.491 km respectively from the reactor cores. The first data set from this experiment was released 
in March 2012 [5] confirming that ^13 was indeed non-zero at 4.9(7 C.L.. We consider in our analysis 
simulated data with 3 years of full run for the RENO experiment and include an uncorrelated 
systematic uncertainty of 0.5% which is the projected benchmark systematic uncertainty for RENO 

• Daya Bay: The Daya Bay reactor experiment [50] observes antineutrinos from the Daya Bay and 
Ling Ao I and Ling Ao II reactors. Each of them have two reactor cores with a total combined 
power of 17.4 GWth- This experiment when fully constructed will have 4 near detectors each with 
20 ton fiducial mass and 4 far detectors also with fiducial mass 20 ton each. The far detectors will 
be at a distance of 1.985, 1.613 and 1.618 km respectively from Daya Bay, Ling Ao and Ling Ao 
II reactors. The distance of the near detectors for each of the reactor cores is more complicated, 
and can be found in [50]. The experiment has been running with 6 detectors and so far produced 
outstanding results [4, 6j. We analyse simulated data corresponding to 3 years of full run for the 
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Figure 1: The efficiencies and resolutions for muons in ICAL@INO as a function of the true muon energy. 
The top left panel shows the reconstruction efficiency of muons, the top right panel shows the charge 
identification efficiency. The bottom left panel shows the zenith angle resolution in cos Q while the 
bottom right panel shows the energy resolution of the muons. The lines in three colors correspond to 
three different benchmark zenith angles for the muons. 



Daya Bay experiment and consider an uncorrelated systematic uncertainty of 0.18% which is the 
projected systematic uncertainty for Daya Bay. 

• T2K: In the T2K experiment [51j a 2.5° off-axis neutrino beam is sent from the J-PARC accelerator 
facility at Tokai to the Super-Kamiokande detector at Kamioka at a distance of 295 km. The beam 
power used for the simulation is taken to be 0.75 MW with 5 years of neutrino running. The fiducial 
mass of Super-Kamiokande is 22.5 kton and there is a near detector ND280 at a distance of 280 m 
from the beam target. 

• NOi^A: The NOz^A experiment |52j will shoot a 3.3° off-axis (anti) neutrino beam from NuMI at 
Fermilab to the 15 kton Totally Active Scintillator Detector (TASD) located in Northern Minnesota 
at a distance of 810 km. The near detector at Fermilab is a 200 ton detector similar to the far 
detector. The beam power used is 0.7 MW with 3 years of running in the neutrino and 3 in the 
antineutrino mode. 
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Figure 2: The number of fi^ events for 10 years of running of ICAL@INO. The left panel shows the 
events in zenith angle bins for muon energy range 3 — 4 GeV. The right panel shows the energy event 
spectrum in the zenith bin —0.6 to —0.5. The long-dashed (red) lines are the distribution of fi~ events for 
Nuance events on which only the oscillations are imposed. The short-dashed (blue) lines give the zenith 
angle distribution when we fold the energy and angle resolution on the earlier event spectrum. The solid 
(green) lines are obtained when in addition to the resolutions we also fold in the efficiencies for muons in 
ICAL@INO. 

4 Simulated Events in ICAL@INO 

For the atmospheric neutrino analysis presented in this paper, we have generated the unoscillated and 
events using the Nuance event generator |32j adapted for ICAL@INC[1. In order to reduce the Monte 
Carlo fluctuations in the event sample, we generate events corresponding to 50 x 1000 kton- years exposure, 
which corresponds to 1000 years of running of ICAL@INO. This event sample is finally normalized to a 
realistic number of years of running of ICAL@INO, when we statistically analyze the data in the end. 
Since it takes fairly long to run the Nuance code to generate such a large event sample, running it over 
and over again for each set of oscillation parameter is practically impossible. Therefore, we run the event 
generator only once for no oscillations and thereafter impose the reweighting algorithm to generate the 
event sample for any set of oscillation parameters. This reweighting algorithm works according to the 
following prescription: 

Every event given by the event generator is characterized by a certain muon energy and muon 
zenith angle, as well as a certain neutrino energy and neutrino zenith angle. For this neutrino energy and 
neutrino zenith angle, the probabilities P^^u^ and Pufj^ue are calculated numerically for any given set of 
oscillation parameters. A random number R between 0-1 is generated. If i? < P^^ue it is classified as a 
z^e event. U R> {Py^u, +Pu^u^), then we classify the event as a Vr event. If P^^^ue < R< (Pu^f^ + Pu^uJ, 

■^We use the Honda et. al atmospheric neutrino fluxes calculated for Kamioka [53]. While the atmospheric neutrino fluxes 
for Theni have been made available recently [54], they are yet to be fully implemented in the ICAL@INO event generator. 
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then it means that this event has come from an atmospheric which has survived as a and is hence 
selected as muon neutrino event. Of course only the muon events are relevant for us, while the others are 
essentially discarded (in this work). Since we do this for a statistically large event sample, we get a i/^j 
"survived" event spectrum that follows the survival probability to a high precision. One could also get 
muon events in the detector from oscillation of atmospheric fg into v^. To find these events we generate 
events from Nuance using atmospheric Ve fluxes but charged current interactions in the detector, with 
the oscillation probability part of the code switched off. In order to get the oscillated muon event from 
in this sample we take a random number S and use a similar procedure for classifying the events. That 
is, if S" < Pueu^^ then the event is taken as an "oscillated" event, which is the only part relevant for us 
in this work. The net number of /i~ events are obtained by adding the "survived" and the "oscillated" 

events. The same exercise is performed for generating the events in the detector. We have checked 
that the final event spectrum we obtain with our method agrees remarkably well with the event spectrum 
obtained by passing the oscillated fluxes through the event generator. 

The data sample after incorporating oscillations is then distributed in very fine muon energy and 
zenith angle bins. This raw binned data is then folded with detector efficiencies and resolution functions 
to simulate the reconstructed muon events in ICAL. In our work we have used (i) the muon reconstruction 
efficiency, (ii) the muon charge identification efficiency, (iii) the muon zenith angle resolution and (iv) the 
muon energy resolution, which have been obtained by the INO collaboration [33]. The energy and zenith 
angle resolutions for the muons are in the form of a two-dimensional table. This means that the muon 
energy resolution is a function of both the muon energy as well as the muon zenith angle. Likewise the 
muon zenith angle resolution depends on both the muon energy as well as muon zenith angle. The muon 
detection efliciencies are also dependent on both the energy and zenith angle of the muon. Similarly 
the charge identification efficiency depends on both the energy and angle of the muon. In Fig. [T] we 
show a snapshot of the efficiencies and resolutions for the ^~ events in ICAL@INO@. They are shown as a 
function of the muon energy for three specific muon zenith angle bins of width 0.1 at cos B = —0.85, —0.55 
and —0.35. The top-left panel shows the reconstruction efficiency of muons as a function of the muon 
energy. The top-right panel shows the charge identification efficiency. The bottom-left panel gives the 
muon zenith angle resolution, while the bottom-right panel shows the muon energy resolution. One can 
notice that there is a rather strong dependence of all the four quantities on the muon energy as well zenith 
angle. The reconstruction efficiency, charge identification efficiency and the zenith angle resolution are 
seen to improve with muon energy. The muon energy resolution on the other hand shows a more complex 
behavior. While for energies 1 — 5 GeV the ge/E is seen to decrease as energy is increased, thereafter 
it increases with energy. In the figure the detector performance is seen to be best for cos O = —0.85 
and worst for cos = —0.35. This is related to the geometry of the ICAL detector wherein there are 
horizontal slabs of iron and RPCs. As a result, the more horizontal muons travel longer in iron and 
hit a lesser number of RPCs. Therefore, the detector performs worse for more horizontal bins. In fact, 
for zenith bins —0.1 ^ cosO ^ 0.1 it becomes extremely difficult to reconstruct the muon tracks and 
hence for these range of zenith angles the reconstruction efficiency is effectively zero. The efficiencies and 
resolutions for and events in ICAL@INO have been obtained separately from simulations and are 
found to be similar [33]. We use the separate and /i"^ efficiencies and resolutions in our simulations 
and results. 

The detector response to muons used in this work comes from the first set of simulations done with the ICAL code. 
These simulations are on-going and are being improved. Hence the values of the resolution functions and efficiencies are 
likely to evolve along with the simulations. More details on this will appear shortly in a separate paper on the detector 
response to muons from the INO collaboration on the ICAL Geant based simulations [33| . 
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The number of /x events in the i bin after implementing the efficiencies and resolutions are given 

C ^max PCOS ©max 

dE / 

^min COS ©min 



N'^^'{fi-)=NY, r^^^dE r"°™dcose Re Re (£iCin,{^l-) + Si{l - Ci)n,{^,+) \ , (5) 



where J\f is the normalization required for a specific exposure in ICAL@INO, and ni{fi~) and ni{fj,~^) are 
the number of and /i"*" events in the i^^ bin respectively, obtained by folding the raw events from Nu- 
ance with the three-generation oscillation probabilities using the reweighting algorithm and subsequently 
binning the data as described earlier. The summation is over i where i scans all true energy and angle 
bins. In Eq. ([5]) Et and cos Qt are the true (kinetic) energy and true zenith angle of the muon, while 
E and cos Q are the corresponding (kinetic) energy and zenith angle reconstructed from the observation 
of the muon track in the detector. In Eq. 8i and fj are the reconstruction efficiencies of /x~ and 
^'^ respectively, while Ci and Ci are the corresponding charge identification efficiencies for /i^ and 
respectively, in the i*^ bin. Both the reconstruction efficiencies as well as charge identification efficiencies 
are functions of the true muon energy Erp and true muon zenith angle cos 0^. In the energy range above 
1 GeV, the resolution functions Re and Rq are seen to be Gaussian from ICAL simulations ^33j and are 
given by 

.[E^-Ef 

4 



Re = Ne exp ( ^ ) ' (6) 



/ - (cos ©T - cos 6)^ \ 
Re = Ne exp (^^ ^ ^ j , (7) 

respectively. A'^^; and Nq are the normalization constants which we numerically calculate and use in our 
code. The resolution functions aE and ae are obtained from ICAL simulations [33] and depend on both 
Et and cos By. An expression similar to Eq. ([5]) can we written for the /x"*" events N^^{ii^). 

In Fig. [5] we show the ^~ event distribution expected in ICAL@INO with 10 years of exposure. In 
the left panel the events are shown for the energy bin 3 — 4 GeV and in cos bins of width 0.1. The right 
panel shows events in the zenith range cos = —0.6 to —0.5 and in energy bins of width 1 GeV. The red 
long-dashed lines shows the Nuance events obtained after including the effect of oscillations according to 
the reweighting algorithm described above. The oscillation parameters used are given in Table [H with 
sin^ ^23 = 0.5, sin^ 2^13 = 0.1 and assuming normal hierarchy. The blue short-dashed lines in the figure 
are obtained once the oscillated events (red long-dashed lines) in ICAL@INO are folded with the energy 
and zenith angle resolutions. A comparison of the red long-dashed and blue short-dashed lines in the 
right panel of the figure reveals that the effect of the energy resolution is to flatten the shape of the 
energy spectrum. We notice that the blue short-dashed line falls below the red long-dashed line for lower 
energy bins, while the trend is reversed for the the higher energy bins. On the other hand, the impact 
of the angle resolution is seen to be negligible for most of the zenith angle bins, as can be seen from 
from left panel of the figure. The reason for these features can be found in the size of the la width of 
the resolution functions shown in Fig. [TJ The energy resolution is seen to be (Je — 0.15E' and hence for 
E between 1 — 11 GeV we do expect some spill-over between bins leading to a smearing of the energy 
spectrum. On the other hand, the bottom right panel of Fig. [T] reveals that the cos resolution for most 
of the zenith angle bins are seen to be better than cos ~ 0.01 — 0.02, while the zenith bins that we 
have used are AcosO = 0.1 in width. This is why the smearing due to angle resolution is essentially 
inconsequential for this choice of zenith angle binning. Note that despite the extremely good detector 
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Figure 3: ICAL@INO /i event spectrum in zenith angle bins for four specific energy bins of E = 4 — 5 
GeV (top left panel), E = 5-6 GeV (top right panel), E = 6-7 GeV (bottom left panel), and S = 7-8 
GeV (bottom right panel) and for 10 years exposure. The solid blue lines correspond to N^^{^~) for 
normal hierarchy while red-dashed lines are for inverted hierarchy. 

performance for angle reconstruction, we still work with such wide angle bins in order to get reasonable 
number of events in each bin. Since the atmospheric neutrino flux falls sharply with energy, the number 
of events per bin becomes very small for very small zenith angle bins in the energy range of 5 — 10 GeV, 
where we expect maximum earth matter effects and hence maximum hierarchy sensitivity. This makes 
the statistical analysis of the data rather challenging. As a compromise we keep bin widths of 0.1 in 
cos 0. 

The green solid lines in Fig. [2] show the realistic events in ICAL@INO where we have taken into 
account oscillations, detector resolutions as well as reconstruction and charge identification efficiencies. 
Meaning, these lines are obtained by imposing the reconstruction and charge identification efficiencies on 
the red short-dashed lines. The left panel shows that once the detector efficiencies are folded, the number 
of events go to almost zero for the horizontal bins. This happens because of difficulty in reconstructing 
the muon tracks along the nearly horizontal directions as discussed before. 

Fig. [3]shows the ICAL@INO /i~ event spectrum in zenith angle bins for four specific muon energy bins 
of ^ = 4 - 5 GeV (top left panel), E = b-Q GeV (top right panel), ^ = 6 - 7 GeV (bottom left panel), 
and E = 7 — 8 GeV (bottom right panel) and for 10 years exposure. The solid blue lines correspond 
to N'j^^{ix~) for the normal hierarchy while red-dashed lines are for the inverted hierarchy. Events were 
generated at the benchmark oscillation point given in Table [H with sin^ 621, = 0.5 and sin^ 2^13 = 0.1. We 
can see that for normal hierarchy there are earth matter effects in the //~ channel leading to suppression 
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of the event spectrum. The extent of the suppression is seen to depend on both energy as well as zenith 
angle bin of the The good energy and angle resolution of the detector is crucial here for fine enough 
binning of the events to extract maximum effect of the earth matter effects. For inverted hierarchy there 
are no earth matter effects in the ^~ events. On the other hand, for /i+ earth matter effects appear for 
inverted hierarchy and are absent for the normal hierarchy. This is why charge separation is so crucial 
for mass hierarchy determination. If the ^~ and /i+ events were to be added one would lose sensitivity 
to earth matter effects and hence to neutrino mass hierarchy. ICAL@INO will have excellent charge 
identification capabilities and zenith angle resolution function as well as good energy resolution. 



5 The Statistical Analysis 

In what follows, we generate the data at the benchmark true values for oscillation parameters given in 
Table [Hand assuming a certain neutrino mass hierarchy. We then fit this simulated data with the wrong 
mass hierarchy to check the statistical significance with which this wrong hierarchy can be disfavored. 
For doing this statistical test we define a for ICAL@INO data as 



Xino 



mm 

{5. 



' 1=1 . 



(8) 



where 



iVf(^")=iVr(^-)|l + ^7r^e, 



(9) 



We have assumed Poissonian distribution for the errors in this definition of x^- The reason is that for 
the higher energy bins ~ 5 — 10 GeV where we expect to see the hierarchy sensitivity, the number of 
events fall sharply (cf. Fig. [3|) and for small exposure times these bins could have very few events per 
bin. Since ICAL@INO will have separate data in and we calculate this xfnoil^ ) Xinoit^'^) 
separately for the fj.~ sample and the fi~^ sample respectively and then add the two to get the xfno 



Xino 



Xino 



~^ Xino 



(10) 



In the above equations, Nf^{^^) and N[^{fi^) are the observed number of /i" and /u^ events respectively 
in the i^^ bin and N[''^[^~) and N-^^{fi'^) are the corresponding theoretically predicted event spectrum 
given by Eq. ([5]) . This predicted event spectrum could shift due to the systematic uncertainties and this 
shifted spectrum Nj^^ is given by Eq. ([9]). In the above irj is the j*'* systematic uncertainty in the i^^ bin 
and is the pull variable corresponding to the uncertainty tt-'. The Xino minimized over the full set 
of pull variables {^j}- The index i runs from 1 to the total number of event bins given by Ni,in- In our 
analysis we have considered the muon energy range 1 GeV to 11 GeV with 10 bins of bin size 1 GeV. 
The zenith angle range in cos0 is taken from —1 to +1 with 20 bins of bin size 0.1. Therefore, in our 
analysis Nf,in = 200. The index j in Eqs. ([5D and Q runs from 1 to k, where k is the total number of 
systematic uncertainties. We have included the following five systematic uncertainties in our analysis. 
An overall flux normalization error of 20% is taken. A 10% error is taken on the overall normalization 
of the cross-section. A 5% uncertainty on the zenith angle dependence of the fluxes is included. An 
energy dependent "tilt factor" is considered according to the following prescription. The event spectrum 
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Figure 4: Left panel shows the ^Xino wrong hierarchy when normal hierarchy is taken to be true, 

while the right panel shows the corresponding reach when inverted hierarchy is taken as true. The three 
lines are for three different values of sin^ 2^13 (true) = 0.08, 0.1 and 0.12 as shown in the legend box, while 
sin^ 023(true) = 0.5 for all cases. We take only ICAL@INO data into the analysis and in the fit keep all 
oscillation parameters fixed at their benchmark true values. 

is calculated with the predicted atmospheric neutrino fluxes and then with the flux spectrum shifted 
according to 

ME) = ME) (;^) ' ^ ME) (l + 6\n^y (11) 

where Eq = 2 GeV and 5 is the lex systematic error which we have taken as 5%. The difference between 
the predicted events rates for the two cases is then included in the statistical analysis. Finally, an over 
all 5% systematic uncertainty is included. 

For the analysis of the data from the current accelerator and reactor experiments we have used the 
individual Xi foi' each experiment as defined in GLoBES [4 7) . 

6 Mass Hierarchy Sensitivity 

In Fig. S] we show the discovery potential of ICAL@INO alone for the neutrino mass hierarchy, as a 
function of the number of years of running of the experiment. The data is generated for the values of the 
oscillation parameters given in Table [1] and for sin^ ^23(true) =0.5. The three lines correspond to the 
different values of sin^ 2^13 (true) shown in the legend boxes in the figure, which we have chosen around 
the best-fit and 2(t range of the current best-fit. The left-hand panel is for true normal hierarchy while 
the right-hand panel is for true inverted hierarchy. These plots show the sensitivity reach of ICAL@INO 
when all oscillation parameters are kept fixed in the fit at the values at which the data was generated. 
For sin^ 2^13 (true) around the current best-fit of 0.1, we can note from the these plots that with 5 years 
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Figure 5: Same as Fig. H] but here oscillation parameters |Amgg|, sin^ ^23 and sin^ 2^13 are allowed to 
vary freely within their current 3a ranges given in Tabled) 

of ICAL@INO data alone, we will have a 1.8a {1.8a) signal for the wrong hierarchy if normal (inverted) 
hierarchy is true. After 10 years of ICAL@INO data, this will improve to 2.5a {2.5a) signal for the 
wrong hierarchy if normal (inverted) hierarchy is true. The sensitivity obviously increases with the true 
value of sin'^ 2^13 (true). The Ax^ is seen to increase almost linearly with exposure. This is not hard to 
understand as the hierarchy sensitivity comes from the difference in the number of events between normal 
and inverted hierarchies due to earth matter effects. Since this is a small difference, the relevant statistics 
in this measurement is small. As a result the mass hierarchy analysis is statistics dominated and one can 
see from Eq. pOj) that in the statistics dominated regime the Axf^o increases linearly with exposure. 

The hierarchy sensitivity quoted above are for fixed values of the oscillation parameters. This ef- 
fectively means that the values of all oscillation parameters are known with infinite precision. Since 
this is not the case, the sensitivity will go down once we take into account the uncertainty in the value 
of the oscillation parameters. The oscillation parameters which affect the mass hierarchy sensitivity of 
ICAL@INO the most are lAm^gl, sin^ ^23 and sin^20i3. In Fig. [5] we show the mass hierarchy sensitiv- 
ity reach of ICAL@INO with full marginalization over the oscillation parameters |Amgg|, sin^ ^23 and 
sin^ 2^13 , meaning these oscillation parameters are allowed to vary freely in the fit within their current 
3a ranges, and the minimum of the taken from that. The CP phase 6c p does not significantly impact 
the ICAL@INO mass hierarchy sensitivity. This will be discussed in some detail later. Therefore, we 
keep 6cp fixed at in the fit. The parameters Am2i and sin^ 612 also do not affect xfno hence 
are kept fixed at their true values given in Table [TJ From the figure we see that for full marginalization 
within the current 3a allowed range for |Amgjj|, sin'^ ^23 and sin^ 2^13, the sensitivity reach of ICAL@INO 
with 10 (5) years data would drop to 2.2o" {1.6a) for sin^ 023(true) = 0.5 and sin^ 2^13 (true) = 0.1, for 
true normal hierarchy. The impact for the inverted hierarchy case is seen to be more. However, this 
is not a fair way of assessing the sensitivity reach of the experiment since all values of the oscillation 
parameters are not allowed with equal C.L. by the current data. This implies that when they deviate 
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Figure 6: The Ax^ for the wrong hierarchy obtained from a combined analysis of all experiments including 
ICAL@INO as well NOi^A, T2K, Double Chooz, RENO and Daya Bay experiments. The left panel is for 
normal hierarchy taken as true in the data while the right panel is for inverted hierarchy as true. The 
three lines are for three different values of sin^ 20i3(true) = 0.08, 0.1 and 0.12 as shown in the legend box, 
while sin^ ^23 (true) = 0.5 for all cases. In the fit we allow all parameters to vary within their 3a ranges 
as shown in Table [TJ 



from their current best-fit value in the fit, they should pick up a from the data of the experiment (s) 
which constrains them. Therefore, one should do a global fit taking all relevant data into account to find 
the correct estimate of the reach of combined neutrino data to the neutrino mass hierarchy. One way to 
take this into account is by introducing priors on the parameters and adding the additional Xprior the 
fit, analogous to what we had explained in Eq. (j4]) and the related discussion for the priors on the solar 
parameters. Moreover, all oscillation parameters are expected to be measured with much better precision 
by the on-going and up-coming neutrino experiments. In fact, by the time ICAL@INO is operational, all 
of the current accelerator-based and reactor experiments would have completed their scheduled run and 
hence we expect that by then significant improvements in the allowed ranges of the oscillation parame- 
ters would have been made. In particular, we expect improvement in the values of sin'^ 2^13, |Amgg| and 
sin^ 2^23 from the data coming from current accelerator and reactor experiments. A projected combined 
sensitivity analysis of these experiments shows that the la uncertainties on the values of sin^ 2^13, |Am^g| 
and sin^ 26*23 are expected to go down to 0.1%, 2% and 0.65%, respectively |35l- Since marginalization 
over these parameters makes a diff^erence to Xino the wrong hierarchy (cf. Figs. |4]and[5]), better 
measurement on them from the current experiments will therefore improve the mass hierarchy sensitivity 
reach of ICAL@INO. As mentioned before, one could incorporate this information into the analysis by 
including "priors" on these parameters. The sensitivity reach of ICAL@INO with projected priors on 
[Am^gl and sin^ 2^23 keeping other parameters fixed can be found in [55]. In the plot presented in [55] 
la prior of 2% on |Am^g| and 0.65% on sin'^ 2^23 was assumed. Note however, that in that analysis only 
two systematic uncertainties were included in the fit, an overall fiux normalization uncertainty of 20% 
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and an overall cross-section uncertainty of 10%. We will discuss the impact of systematic uncertainties 
again in section [71 In this work we improve the analysis by performing a complete global fit of the atmo- 
spheric neutrino data at ICAL@INO combined with all relevant data which would be available at that 
time, viz., data from the full run of the T2K, NOi^A, Double Chooz, Daya Bay, and RENO experiments. 
The combined sensitivity to the neutrino mass hierarchy as a function of number of years of run of the 
ICAL@INO atmospheric neutrino experiment is shown in Fig. [6l For each set of oscillation parameters, 
the joint from all experiments is given by 

x' = xlo + T.xl (12) 
i 

where J2i xf is the contribution from the accelerator and reactor experiments and i runs over T2K, 
NOi/A, Double Chooz, Daya Bay, and RENO experiments. This joint is computed and marginalized 
over all oscillation parameters. The minimized joint Ax^ is shown in Fig. [6l We reiterate that the 
X— axis in this figure shows the number of years of running of ICAL@INO only, while for all other 
experiments we have considered their complete run as planned in their letter of intent and/or Detailed 
Project Report, as mentioned in section 13. 2i The left panel of the figure shows the sensitivity reach 
if normal hierarchy is true while the right panel shows the reach when the inverted hierarchy is the 
true hierarchy. As in Figs. H] and [5] we generate the data at the values of the oscillation parameters 
given in Tableland with sin^ 023(true) = 0.5 and three different values of sin^ 2^13 (true) = 0.08, 0.1, 
and 0.12. The figure shows that inclusion of the accelerator and reactor data increases the sensitivity 
such that with just 5 years of ICAL@INO data one would have more than 2a evidence for the neutrino 
mass hierarchy even if sin 26'i3(true) = 0.08. For the current best-fit of sin^ 2(9i3(true) = 0.1 we would 
rule out the wrong hierarchy at 2.6a while for larger sin^ 2^13 (true) = 0.12 mass hierarchy could be 
determined with almost 3a C.L.. With 10 years of ICAL@INO data the sensitivity would improve to 
2.6a for sin^ 26*13(1^6) = 0.08, 3. la for sin^ 2(9i3(true) = 0.1 and 3.6a for sin^ 26*13 (true) = 0.12. 

The inclusion of the accelerator and reactor experiments into the analysis improves the sensitivity 
reach to the neutrino mass hierarchy in the following two ways. Firstly, inclusion of these data sets into 
the analysis effectively restricts the allowed ranges of oscillation parameters | Am^g | , sin^ 623 and sin^ 2^13 
such that the statistical significance of the mass hierarchy determination from ICAL@INO alone goes up 
to what we were getting in Fig. [3]for fixed values of the oscillation parameters. In addition, we also get a 
contribution to the mass hierarchy sensitivity from the accelerator and reactor experiments themselves. 

We show in Table [2] the separate contributions from the individual experiments to the statistical sig- 
nificance for the mass hierarchy sensitivity from the global fit. The data is generated for normal hierarchy 
and the benchmark true values of the oscillation parameters given in Tabled] with sin^ ^23 (true) = 0.5 
and sin^ 2^13 (true) = 0.1. The oscillation parameters are allowed to vary freely in the fit and the min- 
imum global x^ selected. The Table [2] shows the individual Ax^ contributions from each experiment 
at this global best-fit point for the oscillation parameters, as well as the combined Ax^- The global 
best-fit for the inverted hierarchy corresponds to Am-21 = 7.5 x 10~^ eV^, sin^ 612 = 0.31, sin^ ^23 = 0.5, 
sin^ 26*13 = 0.1, Amlg = -2.4 x 10"^ eV^ and Scp = 252°. Note that (cf. Eq. ([3])) since Am^i depends 
on the value of 6cp, the change of 6cp in the fit gives lAmg^l = 2.34 x 10"'^ eV^ at the global best-fit 
for the inverted hierarchy (Arrigg < 0). However, the data was generated at Am^g = -1-2.4 x 10"'^ eV^, 
which for 6cp{tTue) = 0° gives lAml^^l = 2.44 x 10"'^ eV^. 

Table [2] shows that it is mainly ICAL@INO and NOz^A which contribute to the Ax^ since reactor 
experiments have no sensitivity to the neutrino mass hierarchy and the baseline for T2K is far too short 
to allow for any significant earth matter effect in the signal. The NOi^A experiment on the other hand has 
a baseline of 810 km and a higher energy neutrino beam. This gives the experiment sizable earth matter 
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Expts 


NOvA 


T2K 


DB 


RENO 


DC 


INO 


ALL 


r.n 


2.59 


0.26 


0.53 


0.12 


0.02 


6.16 


9.68 




2.49 


0.31 


0.63 


0.14 


0.02 


6.17 


9.78 



Table 2: Contribution to the Ax^ towards the wrong mass hierarchy at the global best-fit from the 
individual experiments. The normal hierarchy was taken as true and data was generated at the benchmark 
true values of the oscillation parameters given in Tabled] with sin^ 923{tiue) = 0.5 and sin^ 20i3(true) = 
0.1. The first row (^X^„2 ) shows the individual contributions to the statistical significance when Amgg 

cff 

is used for defining the normal (Am^g > 0) and inverted (Am^g < 0) mass hierarchy, while the lower 
row (Ax^^2 ) gives the corresponding contributions when Aml;^ > is taken as normal hierarchy and 

Amgj^ < as inverted hierarchy. 

effects which in turn brings in sensitivity to the neutrino mass hierarchy. The small Ax^ contribution 
from the reactor experiments comes from the fact that the best-fit value of the oscillation parameters, 
and in particular IAtti-Ij^I which controls the spectral shape of these experiments is slightly different at 
the global best-fit for inverted hierarchy, as discussed above. The T2K experiment on the other hand 
returns a small contribution since the best-fit dcp is different from (5c'p(true) = taken in data. 

The first row of Table [2] shows the Ax^ for the case used in this paper where normal hierarchy is 
defined as Am'^g > and inverted hierarchy is defined as Am'^g < 0. In this case the survival probability 
Pufj^u^ is almost the same for the normal and inverted hierarchies when ^13 = 0. However, this is not 
true for the channels Pu^j^ue relevant for T2K and NOvA and P^^i^^ relevant for the reactor experiments. 
Thus in the best-fit, the value of |Am3;^| and 6cp has to be suitably adjusted so that the oscillation 
probabilities in these oscillation channels are closest to each other in the data and in the fit. This results 
in a small Ax^ from experiments that by themselves have practically no hierarchy sensitivity at all. We 
have indeed checked that the Ax^ for mass hierarchy is zero for all reactor experiments and T2K, when 
the fit is performed using data from only one experiment at a time and all oscillation parameters allowed 
to vary in the fit. Only when one performs a combined analysis does this tension between choice of 
the oscillation frequencies arise, returning a small contribution to the mass hierarchy from the reactor 
experiments and T2K. 

In order to check the impact of our choice of Am^g > (Am^g < 0) as normal (inverted) hierarchy, 
we repeated our global fit by taking Arn^i > as the definition for normal hierarchy and Am'^i < as 
the definition for inverted hierarchy. The result for this case is shown in the second row of Table [2j One 
can see that there is hardly any impact of the definition of mass hierarchy on our final results. 

The results in this section including those shown in Table [2] are for 6cp{ticue)=0. However, the 
accelerator-based long baseline experiments are sensitive to the — >■ Ve oscillation channel. The size of 
this Pu^ue oscillation probability and the resultant sensitivity depends crucially on the CP phase 6cp- 
Therefore, the contribution from NOvA to the statistical significance with which we can determine the 
neutrino mass hierarchy will depend crucially on the value of 6cp{tiue). We will discuss this in some 
detail in section [9l 



18 



pSORMALHIERARCHYl INVERTED HIERARCHY 

I = I I I — — 




ICAL@INO years ICAL@INO years 



Figure 7: The impact of systematic uncertainties on mass hierarchy sensitivity. The red hnes are obtained 
without taking systematic uncertainties in the ICAL@INO analysis, while the green lines are obtained 
when systematic uncertainties are included. Long-dashed lines are for fixed parameters in theory as in 
data, while solid lines are obtained by marginalizing over |Amgg|, sin^ ^23 and sin^ 2^13. 

7 Impact of systematic uncertainties 

The atmospheric neutrino fluxes have large systematic uncertainties. In order to study the impact of these 
systematic uncertainties on the projected reach of ICAL@INO to the neutrino mass hierarchy, we show 
in Fig. [7] the mass hierarchy sensitivity with and without systematic uncertainties in the ICAL@INO 
analysis. The Ax^ is shown as a function of the number of years of exposure of the experiment. The data 
was generated at the benchmark oscillation point. The red lines are obtained without taking systematic 
uncertainties in the ICAL@INO analysis, while the green lines are obtained when systematic uncertainties 
are included. The long-dashed lines are for fixed parameters in theory as in the data, while the solid 
lines are obtained by marginalizing over |Amgfj|, sin'^ 623 and sin^ 2^13. The left panel is for true normal 
hierarchy while the right panel is for true inverted hierarchy. The effect of taking systematic uncertainties 
is to reduce the statistical significance of the analysis. We have checked that of the five systematic 
uncertainties, the uncertainty on overall normalization of the fluxes and the cross-section normalization 
uncertainty have minimal impact on the final results. The reason for that can be understood from the 
fact that the atmospheric neutrinos come from all zenith angles and over a wide range of energies. The 
overall normalization uncertainty is the same for all bins, while the mass hierarchy dependent earth 
matter effects, are important only in certain zenith angle bin and certain range of energies. Therefore, 
the effect of the overall normalization errors get cancelled between different bins. On the other hand, 
the tilt error could be used to modify the energy spectrum of the muons in the fit and the zenith angle 
error allows changes to the zenith angle distribution. Therefore, these errors do not cancel between the 
different bins and can dilute the significance of the data. In particular, we have checked that the effect of 
the zenith angle dependent systematic error on the atmospheric neutrino fluxes has the maximum effect 
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Figure 8: Same as Fig. [Hbut for sin^ 023(true) 
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on the lowering of the Ax^ for the mass hierarchy sensitivity. 

8 Impact of sin^ ^23(true) 

It is well known that the amount of earth matter effects increases with increase in both ^13 and 623- In 
the previous plots, we showed the mass hierarchy sensitivity for different allowed values of sin^ 2^13 (true), 
while sin^ 023 (true) was fixed at maximal mixing. In Figs. [8] and [9] we show the sensitivity to the 
neutrino mass hierarchy as a function of number of years of running of ICAL@INO for different values 
of sin^ 2^13 (true) as well as sin^ 023(true). In Fig. [8] we show the Ax^ corresponding to ICAL@INO 
alone and with oscillation parameters fixed in the fit at their true values. This figure corresponds to 
Fig. [3] of the previous section but now with two other values of sin^ 023(true). In Fig. [9] we give the 
combined sensitivity to mass hierarchy of all accelerator and reactor experiments combined with the data 
of ICAL@INO. We reiterate that the x-axis of the Fig. [9] shows the exposure of ICAL@INO, while for 
all other experiments we have assumed the full run time as discussed in section 13. 2[ The red bands in 
Figs. [8] and [9] correspond to sin^ 023(true) = 0.6 while the green bands are for sin^ 023(true) = 0.4. The 
width of each of the bands is mapped by increasing the value of sin^ 20i3(true) from 0.08, through 0.1, 
and up to 0.12. As seen in the previous subsection, the Ax^ for the wrong mass hierarchy increases with 
sin^ 2^13 (true) for a given value of sin^ 023(true) and ICAL@INO exposure. A comparison of the Ax"^ for 
different values of sin^ 023(true) reveals that the Ax^ also increases with sin^ 023(true). 

From Fig. [9] one infers that for 5(7p(true) = 0, a combined analysis of all relevant experimental 
data including 5 years of ICAL@INO exposure would give the neutrino mass hierarchy from anywhere 
between 1.8a to 3.8cj, depending on the values of sin2 6l23(true) and sin^ 26*13 (true). With 10 years of 
running of ICAL@INO this would improve to 2. la to 4.5(j, depending on what value of sin^ 2^13 (true) 
and sin^ 023(true) have been chosen by mother Nature. Here we have allowed sin'^ 023(true) to vary 
between [0.4 — 0.6] and sin^ 20i3(true) between [0.08 — 0.12]. We next look at the impact of (5cp(true) on 
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Figure 9: Same as Fig. [6] but for sin^ 023(true) = 0.4 (green band) and sin^ 023(true) = 0.6 (red band). 

the prospects of determining the neutrino mass hierarchy. 

9 Impact of 5c p and (5cp(true) 

So far we had taken (^c'p (true) =0 in the data and varied 5c p in the fit only for the long baseline experi- 
ments. For the analysis of the ICAL@INO data, we had kept 5cp fixed to zero in both the data and the 
theory. The reason was that while the for T2K and NOz^A are strongly dependent on 5c p, the mass 
hierarchy for ICAL@INO shows a very mild dependence on it. We show this dependence explicitly 
in Fig. [10] for 10 years exposure in ICAL@INO and compare it with the corresponding dependence of 
NOi/A (see also [31j)|. We generate the data for normal hierarchy and at the benchmark values of the 
oscillation parameters from Table [1] and with sin^ ^23(ti'ue) = 0.5 and sin^ 2^13 (true) = 0.1. In the fit 
with inverted hierarchy, we keep all oscillation parameters fixed, except 5c p which is varied over its full 
range [0 — 2-k\. The corresponding Ax^ is plotted in Fig. [10] as a function of the 5c p in the fit. The 
red long-dashed line shows the 5c p dependence of Ax^ for ICAL@INO, while the blue solid line shows 
the wild fiuctuation of A^^ expected for NOz^A. Amongst the accelerator and reactor experiments we 
show only NOi^A in this figure as the leading contribution to the mass hierarchy comes from this ex- 
periment. We reiterate that at each point we use the data generated at 5c'p(true)=0. The figure shows 
that when we fit the data with inverted hierarchy, the Ax^ for NOz/A changes from more that 26 for 
5c p — 80° to less than 3 for 5c p — 260°. When marginalized over 5c p in the fit, obviously it will return 

^This figure was shown in fST. However, the analysis in [31] was in terms of neutrinos done with some assumed values of 
the detector resolutions and efficiencies. Since we do here the complete analysis of the ICAL@INO projected data in terms 
of the detected muons and with realistic detector resolutions and efficiencies obtained from ICAL simulations, we reproduce 
a similar plot for completeness. 
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Figure 10: The change of A^^ for the wrong hierarchy as a function of the 6cp chosen in the fit. The data 
was generated for normal hierarchy and 6cp(tTue)=0. The other oscihation parameters in both data and 
theory are fixed at their benchmark values given in Table [H The sohd blue hne shows the change in Ax^ 
with 5cp for NOi^A, while the dashed red line shows the corresponding variation of Ax^ for ICAL@INO. 

the lowest value of the Ax^, which in this case would be 2.59 [^l. When marginalized over all oscillation 
paremeters, the sensitivity further reduces to Ax^ = 1.77. The contribution from ICAL@INO on the 
other hand is seen to be almost independent of 6c p- Note that such a figure was also shown in ^55j for 
the ICAL@INO simulations. However, the analysis of the ICAL data has been improved since then and 
more types of systematic uncertainties introduced. This explains the change in the behavior of this plot 
with the dependence of the Ax^ becoming flatter with 5c p- From this study we conclude that one does 
not need to marginalize over 6c p for the ICAL@INO data. However, for the long baseline experiments a 
fine marginalization over this parameter is absolutely crucial, especially for mass hierarchy studies. 

We had seen in Fig. [10] (and as is well known) that the long baseline experiments are very sensitive 
to 6c p- In that figure we were studying the impact of changing 6c p in the fit for a particular 6c p (true) 
in the data. In particular, we had taken 6cp{true)=0. A pertinent question at this point is the following: 
what is the impact of 6c p (true) on the sensitivity of the experiments to the neutrino mass hierarchy? We 
present the answer to this question in Fig. [TT\ where we show the Ax^ for the mass hierarchy sensitivity 
as a function of 6cp{tiue). To obtain these curves, we generate data for normal hierarchy at each value 
of (5cp (true) shown in the j;-axis and then fit this data for inverted hierarchy by marginalizing over 
all oscillation parameters, including 6cp- The data were generated for the benchmark oscillation point 
given in Table [H sin^ 023(true) = 0.5 and sin^ 2^13 (true) = 0.1. The exposure considered is 10 years of 
ICAL@INO and full run for all other experiments. The green solid line in this figure is for only NOz^A, 
the blue short-dashed line is obtained when we combine NOz^A, T2K and all the reactor data, while 

®Note that here the other parameters are fixed and only NOi^A data is being considered in the analysis, while the Ax^s 
shown in Table [5] are for a global fit of all data with all oscillation parameters allowed to vary freely in the fit and the 
combined marginalized over them. 
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Figure 11: Impact of (5cp (true) on the mass hierarchy sensitivity. The sensitivity change of NOz^A due 
to 5c p (true) is shown by the green solid hne. The Ax^ for the wrong mass hierarchy expected from the 
combined data from T2K, NOj^A and the reactor experiments is shown by the short-dashed blue line. 
The global Ax^ for the wrong mass hierarchy from the combined ICAL@INO plus the accelerator and 
reactor data is shown by the red long-dashed line. Data was generated for normal hierarchy and at the 
benchmark oscillation point from Table [1] with sin^ 023 (true) = 0.5 and sin^ 2^13 (true) =0.1 and at each 
value of 5c'p(true) shown in the x-axis. The fit to the wrong inverted hierarchy was fully marginalized 
over all oscillation parameters. The ICAL@INO exposure was taken as 10 years. 

the red long-dashed line is what we get when the ICAL@INO data is also added to the long baseline 
and reactor data. As expected, the ICAL@INO data is almost completely independent of (5cp(true) 
and so is its projected sensitivity to the neutrino mass hierarchy. On the other hand, the reach of the 
NOz^A experiment for determining the neutrino mass hierarchy is seen to be extremely sensitive to the 
value of 5c'p(true). All our plots shown so far on the global mass hierarchy reach were done assuming 
5c'p(true) = 0. We can see from the figure that indeed the statistical significance with which we could 
rule out the inverted mass hierarchy in this case is S.llcr, as discussed before. The Ax^ for 6cp(tTue) = 
for NOz^A is 1.77. However, this quickly falls to almost zero for (5cp (true) ~ [50° — 150°]. Thereafter, it 
rises sharply giving a Ax^ = 8.21 around 5cp(true) ~ 270°, and then finally falls back to Ax^ = 1.77 for 
5c'p(true) = 360°. When T2K and all reactor data are added, there is an improvement to the combined 
sensitivity due to constraint coming from the mismatch between the best-fit for different experiments. 
This is specially relevant in the (5c'p(true) ~ [50° — 150°] range where NOi^A by itself gives no mass 
hierarchy sensitivity. However, once we add the T2K and reactor data to NOi^A data, the Ax^ of this 
combined fit in this region of 5cp (true) increases to ~ 3.5. The reason for this can be understood as 
follows. For the case where 5c'p(true) = 72°, the best-fit for NOi^A alone was dcp = 234°, sin^ 623 = 0.5 
and sin'^ 2^13 = 0.1. For this 6cp{trvLe), T2K data taken alone gave Ax^ ~ with best-fit at 6cp = 198°, 
sin^ ^23 = 0.52 and sin^ 2^13 = 0.08. A combined fit with all accelerator and reactor data gave best-fit 
at 6cp = 198°, sin^ ^23 = 0.48 and sin^ 2^13 = 0.1. This results in a contribution to the mass hierarchy 
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Ax^ = 0.92 from NOj^A, Ax^ = 1-41 from T2K and Ax^ = 1.1 from reactors. That the reactor data 
return a Ax^ contribution to the mass hierarchy sensitivity might appear strange at the outset since the 
combined fit given above has a best-fit sin^ 2^13 = 0.1, and one usually does not expect the reactor data 
to depend on sign of Am'^i, 6cp and sin^ 623. However, note that the reactor data depend on lAm^^l and 
what we use in our fits is Am^g given by Eq. ([3|) which is related to 6cp- Therefore, as discussed earlier 
in section 6, this subtlety regarding the choice of the definition of the neutrino mass hierarchy results 
in a small change in the best-fit IAttt-Ij^I, which in turn results in a small Ax^ contribution to the mass 
hierarchy from the reactor experiments. Finally, addition of the ICAL@INO data raises the Ax"^ by a 
constant amount for all values of 6cp{true). Therefore, depending on (5cp(true) the combined sensitivity 
to the neutrino mass hierarchy could range from 3a (for (5cp (true) = 144°) to 4.2(T (for Jc'p (true) ~ 270°). 
These numbers are for sin^ 023(true) = 0.5 and sin^ 20i3(true) = 0.1 and will improve for larger values of 
these parameters. 

10 Conclusions 

In this paper we looked in detail at the prospects of determining the neutrino mass hierarchy with the data 
collected in the atmospheric neutrino experiment ICAL@INO. The atmospheric muon neutrino events 
before oscillations were simulated using the Nuance based generator developed for ICAL@INO. To reduce 
Monte Carlo fluctuations 1000 years of exposure was used for generating the events. Since it takes very 
long for the generator to produce such a large event sample, we simulated the atmospheric events using 
the generator just once for no oscillations and used a reweighting algorithm to obtain the oscillated event 
sample for any set of oscillation parameters. The oscillated muon event spectrum was then folded with 
the muon reconstruction efficiencies, charge identification efficiencies, energy resolution and the zenith 
angle resolution functions obtained from ICAL simulations (which will appear elsewhere [33]) to obtain 
the reconstructed muon event spectrum in the detector. As ICAL is a magnetized calorimetric detector 
allowing an identification of fi~ and //^ events, it has an edge over rival atmospheric neutrino experiments. 
We defined a x^ function for Poissonian distribution for the errors in the ICAL@INO experiment taking 
into account systematic uncertainties expected in the experiment. 

The data was generated for benchmark true values for the oscillation parameters and a given neutrino 
mass hierarchy and fitted with the wrong hierarchy. We showed the mass hierarchy sensitivity results 
with only ICAL@INO data for the analysis with fixed values of the oscillation parameters in the fit, 
as well as that obtained after marginalization over lAm^gl, sin^ ^23 and sin^ 2^13 in their current 3a 
ranges. We showed these results as a function of the exposure in ICAL@INO. From a comparison of 
the two results, we showed that the mass hierarchy sensitivity with ICAL@INO data deteriorates with 
the uncertainty in the measured value of |Amgg|, sin^ ^23 and sin^ 2^13. These parameters will be rather 
accurately determined by the T2K, NOi^A Double Chooz, RENO and Daya Bay experiments. Since 
INO is expected to start operation after these have finished their full projected run, it is meaningful to 
include their effect in a combined statistical analysis for the neutrino mass hierarchy. In order to take 
that into account, we simulated the data for these experiments using GLoBES with the experimental 
specifications mentioned in their respective Letter Of Intent and/or Detailed Project Report. The results 
on mass hierarchy sensitivity from the combined analysis of data from ICAL@INO along with that from 
T2K, NOi^A, Double Chooz, RENO and Daya Bay was shown for benchmark values of the oscillation 
parameters and full marginalization over all oscillation parameters in the fit for the wrong mass hierarchy. 
We showed that marginalization over 6cp is practically unessential for the ICALCDINO data. However, for 
the accelerator data it is absolutely crucial to marginalize over 5c p due to the very strong dependence of 
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the hierarchy sensitivity on this parameter in these experiments. We then generated the data at all values 
of (5cp (true) and showed that the mass hierarchy sensitivity of ICAL@INO was independent of (^c'p(true), 
however, the sensitivity of the combined NOi^A, T2K and the reactor experiments depends very strongly 
on what 5c'p(true) has been chosen by Nature. For sin^ 023(ti'ue) = 0.5 and sin^ 2^13 (true) = 0.1 the 
combined data of 10 years exposure in ICAL@INO along with T2K, NOz^A and reactor experiments 
could rule out the wrong hierarchy with a statistical significance of 3a to 4.2(T, depending on the chosen 
value of 6cp{tTue). We also studied the effect of sin^ 2^13 (true) and sin^ 023(true) on the reach of these 
combined projected data sets to determining the neutrino mass hierarchy. For Jc'p(true) = 0, we showed 
that the statistical significance with which the wrong hierarchy could be ruled out by the global data 
set comprising of 10 years exposure in ICAL@INO along with T2K, NOz^A and reactor experiments, 
could be anywhere between 2.13(7 to 4.5ct depending on sin^ ^23(true) and sin^20i3, where we allowed 
sin'^ 023(ti"ue) to vary between [0.4 — 0.6] and sin^ 20i3(true) between [0.08 — 0.12]. For the most favorable 
choice of (5c'p(true) ~ 270° the sensitivity could go up to greater than 5a with 10 years of ICAL@INO 
combined with data from T2K, NOi^A and reactor experiments. 
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